module timestep
use technical
!use grid
use mod_2dflu
use equations
implicit none
public 
save
! -----------------------------------------------------------------------------------------------------!
! Timesteps 2-d arrays. Solves the equation
! (d/dt) f = df   where both f and df are two dimensional arrays.
! Actually both f and df has a third dimension which is number of variables
! we use. This is called "tvar".  For example for two dimensional hydrodynamics
!  this is actually one (only vorticity). If we include passive scalar this would
! be two. This subroutine does not care what those variables are. It is just
! given an f and a df and then timesteps. 
! -----------------------------------------------------------------------------------------------------!
character (len=labellen) :: tscheme
double precision, allocatable, dimension(:,:,:) :: fzero,fk1,fk2,fk3
logical :: lallocated_fzero,lallocated_fk1,lallocated_fk2,lallocated_fk3
namelist /timestep_param/ & 
      tscheme,delta
contains
!**************************************************************!
subroutine initialize_timestep
!
! Does initial setup depending on what kind of time stepping scheme
! we are going to use. 
!
!
! First set all the lallocated variables to zero just to make sure that we start with
! right conditions.
!
  lallocated_fzero=.false.;lallocated_fk1=.false.
  lallocated_fk2=.false.;lallocated_fk3=.false.
!
! Then choose the appropriate time step scheme and initialize accordingly.
!
  select case (tscheme) 
  case ('adbsh2')
    call announce('initialize_timestep:','using Adams-Bashforth exponential (2nd order) scheme')
    call initialize_adbsh2
  case ('rk2')
    call announce('initialize_timestep:','using Runge-Kutta (2nd order) scheme')
    call initialize_rnkt2 
  case ('rk4')
    call announce('initialize_timestep:','using Runge-Kutta (4th order) scheme')
    call initialize_rnkt4
  case default
    call fatal_error('initialize_timestep:','time stepping scheme not found')
  endselect
!
endsubroutine initialize_timestep
!**************************************************************!
subroutine initialize_adbsh2
!
! Initializes the Adams-Bashforth (2nd order) time stepping scheme
!
  lexponential_integrator=.true.
endsubroutine initialize_adbsh2
!**************************************************************!
subroutine initialize_rnkt2
!
! As this is a second order scheme there is need to have intermediate storage
! which must be allocated now. 
!
  lexponential_integrator=.false.
  allocate(fzero(n1+2,n2,tvar))
  lallocated_fzero=.true.
!
endsubroutine initialize_rnkt2
!**************************************************************!
subroutine initialize_rnkt4
!
! As this is a second order scheme there is need to have intermediate storage
! which must be allocated now. 
!
  lexponential_integrator=.false.
  allocate(fzero(n1+2,n2,tvar))
  lallocated_fzero=.true.
  allocate(fk1(n1+2,n2,tvar))
  lallocated_fk1=.true.
  allocate(fk2(n1+2,n2,tvar))
  lallocated_fk2=.true.
  allocate(fk3(n1+2,n2,tvar))
  lallocated_fk3=.true.
!
endsubroutine initialize_rnkt4
!**************************************************************!
subroutine read_timestep_param (fpointer)
  integer :: fpointer,ierr
  read(fpointer,NML=timestep_param, IOSTAT=ierr)
  if(ierr/=0) call fatal_error('timstep','problem with reading namelist')
endsubroutine read_timestep_param
!**************************************************************!
subroutine ddt_equ (f,df)
!
! Time stepping the actual equations. 
! Calls different time-steppers depending on which tscheme is being used.
!
  double precision, dimension(n1+2,n2,tvar) :: f
  double precision,dimension(n1+2,n2,naux) :: df
!
  select case (tscheme) 
  case ('adbsh2')
    call adbsh2(f,df)
  case ('rk2')
    call rnkt2(f,df) 
  case ('rk4')
    call rnkt4(f,df) 
  case default
    call fatal_error('ddt_equ:','time stepping scheme not found')
  endselect
!
endsubroutine ddt_equ
!**************************************************************!
subroutine adbsh2 (f,df)
!
  double precision, dimension(n1+2,n2,tvar) :: f
  double precision,dimension(n1+2,n2,naux) :: df
  integer :: istep
  double precision :: tzero
!
endsubroutine adbsh2 
!**************************************************************!
subroutine rnkt2 (f,df)
!
!Numerical Recipes for Fortran Chapter 16, page 704 
! fzero stores the initial value of f.  
!
  double precision, dimension(n1+2,n2,tvar) :: f
  double precision,dimension(n1+2,n2,naux) :: df
  integer :: istep
  double precision :: tzero
!
  fzero(:,:,:) = f(:,:,:)
  tzero=time
  istep=1
  call get_df (istep)
  f(:,:,:) = fzero(:,:,:) + (delta/2.)*df(:,:,1:tvar)
  time=tzero+delta/2.
  istep=2
  call get_df (istep)
  f(:,:,:) = fzero(:,:,:)+delta*df(:,:,1:tvar)
  time=tzero+delta
!
endsubroutine rnkt2
!**************************************************************!
subroutine rnkt4 (f,df)
!
!Numerical Recipes for Fortran Chapter 16, page 705 
!
  double precision, dimension(n1+2,n2,tvar) :: f
  double precision,dimension(n1+2,n2,naux) :: df
  integer :: istep
  double precision :: tzero
!
  fzero(:,:,:) = f(:,:,:)
  tzero=time
  istep=1
  call get_df(istep)
  fk1(:,:,:) = df(:,:,1:tvar)
  f(:,:,:) = fzero(:,:,:) + (delta/2.)*df(:,:,1:tvar)
  time=tzero+delta/2
  istep=2
  call get_df(istep)
  fk2(:,:,:) = df(:,:,1:tvar)
  f(:,:,:) = fzero(:,:,:) + (delta/2.)*df(:,:,1:tvar)
  time=tzero+delta/2
  istep=3
  call get_df(istep)
  fk3(:,:,:) = df(:,:,1:tvar)
  f(:,:,:) = fzero(:,:,:) + delta*df(:,:,1:tvar)
  time=tzero+delta
  istep=4
  call get_df(istep)
  f(:,:,:) = fzero(:,:,:) + (delta/6.)*fk1(:,:,:) + (delta/3.)*fk2(:,:,:) + & 
                                        (delta/3.)*fk3(:,:,:) + (delta/6.)*df(:,:,1:tvar)
  time=tzero+delta
!
endsubroutine rnkt4
!**************************************************************!
subroutine cleanup_timestep
  if(lallocated_fzero) deallocate(fzero)
  if(lallocated_fk1) deallocate(fk3)
  if(lallocated_fk2) deallocate(fk2)
  if(lallocated_fk3) deallocate(fk1)
endsubroutine cleanup_timestep
!**************************************************************!
endmodule timestep
